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The determination of the two-body density functional from its one-body density 
is achieved for Moshinsky's harmonium model, using a phase-space formulation, 
thereby resolving its phase dilemma. The corresponding sign rules can equivalently 
be obtained by minimizing the ground-state energy. 

I. INTRODUCTION 

The harmonium model originally proposed by Moshinsky 1 has earned its spurs as a simple 
analogue to a two-electron atom, helpful to illustrate the main ideas of reduced density 
matrix and correlation energy theory in an exactly solvable context. The model consists 
of two spin- 1 fermions trapped in a harmonic potential and repelling each other with a 
Hooke's law force, as well. Chapter 2 of the book by Davidson 2 describes its ground state 
in the standard wave function formalism, as well as the reduced density matrix and the pair 
distribution, exhibiting correlation. 

A one- dimensional version of the model was put to work by Neal 3 in the hope of finding an 
exact universal density functional of the Hohenberg-Kohn-Sham type. 4 This proves illusory; 
but the computations by Schindlmayr in his very pedagogical rejoinder 5 make it clear that 
Neal's harmonium scheme supports successful approximations for confining potentials. More 
recently, the harmonium model has proved its worth in suggesting approximate general forms 
of 1-density matrices 6 and Ansatze for correlation energy density. 7 

The advantage of the harmonium model is that the required computations can be ana- 
lytically performed. However, despite this solvable character, several pertinent functionals 
have not been exploited so far. It is well known that possession of the 1-body matrix p\ for 
an iV-electronic system does not allow effective inference of the corresponding 2-body ma- 
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trix p 2 , which would trivialize the energy functional in quantum chemistry. It is natural to 
diagonalize p\ and seek to expand p 2 in terms of eigenfunctions fj of p 1 ( "natural orbitals" ) 
and its eigenvalues < rij < 1 ("occupation numbers"), with J2j n j = 1- O ver the years, 
starting with the work by Miiller, approximate functionals based on this spectral analysis 
of p\ have been suggested and tried with various results. 

Two-electron atoms constitute the exception to our ignorance. In this article we focus 
on the exact Shull-Lowdin-Kutzelnigg (SLK) functional for the ground state of such atoms 
in terms of natural orbitals. 2 ' 8 ' 9 Work by those authors in the late fifties and early sixties 
established that, for a reduced 1-density of the kind 

Pl ( X ,x') = iflyiv +uu>)pi(ry) = ^(titi< +u<) E^WJV), 

j 

the corresponding 2-density matrix is given by the form 



P2(X\, X 2 , x[, X 2 ) — -(Tli 2 — ^lt 2 ) (tl'4-2' — -^l'lV) 
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Alas, the SLK recipe, although exact, is underdetermined: of the Cj we only know that 
\cj\ 2 = rij. This is a "phase dilemma" of density functional theory. We work here only with 
states described by real wavefunctions — still leaving us with an infinite number of signs to 
account for. 

Notwithstanding its venerable age, formula (1) apparently has never been verified exactly. 
A theorem without an example is a sorry thing. Of course, numerical computations tend to 
confirm the SLK theorem; but one should not forget that they tell us about the approxima- 
tions (nearly always from a Hartree-Fock starting point), rather than the true solution. We 
verify the SLK method for harmonium in full detail, including the energy functional, in the 
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following three sections. Along the way, we solve the sign conundrum for the model. Our 
methods are elementary, asking familiarity with little more than orthogonal polynomials at 
the level of Lebedev 10 or Andrews et al. u 

Nevertheless, within the standard formalism it is not at all obvious how to go about the 
problem. We manage to sidestep difficulties by working with the Wigner quasiprobability 
on phase space instead. A recent quantum phase space view of harmonium, dealing with 
other matters, is given by Dahl. 12 

In the concluding Section 5 we very briefly discuss the new perspectives on correlation 
energy and approximate functionals for p 2 revealed by the treatment in this paper. 

We follow Davidson's notation 2 as much as feasible. A good review on the Miiller func- 
tional is found in Ref. 13. One may consult Refs. 14, 15 for popular variations on it. 



II. THE SETUP 

The Hamiltonian for harmonium in Hartree units is 

tj Pi . V2 i ^/ 2 i 2\ ^2 /ri\ 

H = y + y + ^ l + r 2 )- -r 12 . (2) 
Introduce extracule and intracule coordinates, respectively given by 

R = (ri + r 2 )/v / 2, r = (n - r 2 )/V2, 
with conjugate momenta 

P=( Pl +p 2 )/V2, p=( Pl -p 2 )/V2. 

Therefore 

u u u P 2 u 2 R 2 p 1 /xV 

H = H R + H r = 1 h — + - — . 

2 2 2 2 

As advertised, our notation is that of Davidson 2 except that our 5 is equal to twice his a, 
and we introduce the frequencies u = y/k and /J, = \Jk — 5; assume < 5 < k for both 
particles to remain in the potential well. 

Since the spin factors are known, we concentrate on the spinless part of the quantum 
states henceforth. The spinless Wigner quasiprobability (normalized to one) corresponding 
to a (real) 2-particle wave function \1> is given by 

iVri, r 2 ; Pl ,p 2 ) = ^J p 2 (n - Zl , r 2 - z 2 ; n + Zl , r 2 + z 2 ) e «(pi-*i+w-«) £ Zx Sz 2 , (3) 
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with p 2 {ri,r 2 ,r' l ,r' 2 ) = ^(r 1 ,r 2 )^(r' 1 ,r' 2 ). The definition extends to transition matrices 
| $)($'| also: 

*W(ri, r 2 ; Pl ,p 2 ) = J $(n - Zi, r 2 - z 2 )$'(n + *i, r 2 + z 2 ) e^ 1 '* 1 ^ 2 ^ rf 3 ^ d 3 z 2 . 

Fourier analysis easily provides the inverse formula to this, that we do not bother to write. 

By use of (3) and the ground state wave function for harmonium, one can obtain the 
corresponding Wigner function, which factorizes into extracule and intracule parts: 



p ( ^ 1 f 2H A ( 

P gs (r-i,r 2 ;pi,p 2 ) = ^expl —J expl - 



2H r 



(4) 

This is reached more efficiently and elegantly by the methods of phase space quantum 
mechanics. 16 One can now obtain p 2 , given by the inverse formula of (3). The pairs density 
P2( r i, r 2, r i, 7*2) is recovered by integration over the momenta. 

The reduced 1-body phase space (spinless) quasidensity for the ground state d gs is ob- 
tained, as in the standard formalism, by integrating out one set of variables, 

d ( r -n)-—f 4UfI V /2 r -2r 2 ^/(^+M) r -2p 2 /^+ M ) (r) 

We leave it as an exercise for the reader to recover Eq. (2-68) of Ref. 2 for pi(ri,r[) from 
this. The marginals of d gs give the electronic density and momentum density. 

It should be recognized that, while P gs is a pure state, mathematically d gs describes a 
mixed state. For Gaussians on phase space, such as P gs and d gs too, there are simple rules 
to determine whether they represent a pure state, 17 a mixed state, 18 or neither. Writing 
q = (n,r 2 ), 7T = (Pi,p 2 ), u = (qr,7r), we find P gs {u) = 7i- 6 e - u - Fu = ^e-i-M-^-A-^ 
where, amusingly, 



2 



00 — n u + p 




We see that the matrix F corresponding to formula (4) is symmetric and symplectic, and 
therefore represents a pure state. This is not the case for d gs . Thus recovering P gs from 
knowledge of d gs alone is akin to putting Humpty Dumpty together again! 

III. COMPUTATION OF THE 2-BODY QUASIDENSITY 

Since all the relevant quantities factorize, in this section we work in one dimension (instead 
of three) for notational simplicity. The real quadratic form in the exponent of d gs must be 
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symplectically congruent to a diagonal one. 18 We perform the transformation 

(Q,P) '■= ((w/i) l74 r, (o>//) _1//4 p); or, in shorthand, U = Su, 

where now u = (r,p). Here S being symplectic just means having determinant 1, which is 
evidently the case. Introducing as well the parameter A := (u+fj), the 1-quasidensity 

takes the simple form 

d gs (u(U)) = -e~ xu2 . 

71 

We may also write A =: tanh(/3/2), so that 

I 3 = log i \ = 2 lo § ~r= — V' and Smh (/V 2 ) = R ro = — • 

1 - A ytuj-^fjl Vl - A 2 to-pL 

From the series formula, valid for \t\ < 1, 



CO 

x/2 f n = 1 r -x(l+t)/2(l-t) 



n=0 

taking t = -(1 - A)/(l + A) = -e _/3 and x = 2U 2 , it follows that 



^sinh|^(-irL r (2f/ 2 )e- c/2 e- 



7T 

r=0 

We recognize the basis of Wigner eigenfunctions on phase space standing for the oscillator 
states: 16 

frr(V) = -(-l) r L r (2U 2 )e- u2 . 

71 

Note the normalization / f 2 r (Q, P) dU = (2tt)~ 1 . Consequently, we realize that d gs is in thin 
disguise a Gibbs state, 18 with inverse temperature j3: 

d gs (u) = d^S-'U) = 2 sinh I J2 e-^ +1 ^ 2 f rr (U). (7) 

r=0 

Thus we have identified the natural orbitals in the U variables. Their occupation numbers 
are 

n r = 2 sinh t e-W* = ^ ( ^Z^B) ^ = ^ ( V^WR ) * (8) 

Notice that no = 1 — e" 13 = where Z is the partition function for the system; also 

J2r n r = (1 ~~ e " /? ) Sr e_r/? = 1- These n r have nice square roots: 



in,, = 
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We prepare now to test the SLK functional. On phase space, formula (1) is replaced by 

j oo 

^2SLk(«1,«2; SPin) = - (tl^2 - I1T2) (tl42' - ll'tv) ^2 Cr ° s frs(Ui)f rs {u 2 ). (9) 

r,s=0 

The f rs are Wigner eigentransitions, the functions on phase space corresponding to matrix 
transitions between oscillator states. They are well known. 16 For r > s, abusing notation, 

/„(«) := -(-l) s \/4 {2U 2 t~ s ^ 2 e~^- s ^ L r ~ s {2U 2 ) e~ u \ 

where d := arctan(P/<3). Then f sr is the complex conjugate of f rs . In (9) we proceed to 
sum over each sub diagonal, where r — s = I > 0: 



^2 ^ n r1T'sfrs{Ul)frs{u2) 



r—s=l 

oo 



^ e^' 2 {2U 1 U 2 ) l e- a ^ + ^ e'^- ]T — ^— L l s (2U 2 )L l s (2U 2 ) 

s=0 ^ ' 



- 1 c -(u?+ul)/x c -Wi+02) L ( 2U ^ \ 
"tt 2 l \smh(P/2))' 

where Ii denotes the modified Bessel function, on use of another series formula: 10 

V n! L a (x)L a (v) t n = e -(**0t/(i-t) / ( 

^An + a)\ n[ j n[V) 1-t 6 la \i- t )- 

n=0 v ' N ' 

Similarly for r — s = —I < 0, we get the same result replaced by its complex conjugate. 
Borrowing finally the generating function identity for Bessel functions, 



I (z)+2j2l l (z) cos(W) = e z 



„Z COS 

1=1 

where, by taking 9 = ■d 1 + # 2 + tt, one obtains for the total sum: 



2 p -[(U?+Ul)/\+2UiU 2 csch(/3/2) cos(tfi+tf 2 )] 



7T e 



- 2 p-5[('??+9l)(^+A')+(p?+p|)(^" 1 +^ 1 )] p -9i92(w-M) p PlP2(^" 1 ^- 1 ) 



which in view of (4) is the correct result. Clearly the choice 9 = $i + d 2 + 7r amounts to the 
alternating sign rule for the functional: 

c r = (— l) r ^/rv for r — 0,1,2,... 
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In the end, for -P 2 slk(wi, u i\ spin) we obtain: 

2(tli 2 -Ilt 2 )(tl4 2 ' -Il't 2 ') ^(-) nr+n3 V / ^/r S (Ml)/r S (M 2 ). 

r,s=0 

As far as we know, this is the first time that the solution to the sign dilemma has been 
exhibited for any model. No big deal, a devil's advocate might say, since P gs was known 
beforehand. But, in point of fact, the correct choice of signs may instead be chosen by 
optimization of the energy functional; so it can be found without prior knowledge of the 
system's ground state. Our next step is to confirm this claim. 



IV. THE ENERGY FUNCTIONAL 

We still work in dimension one. The energy E gs of the ground state is of course u/2 + /x/2. 
This contains one-body contributions Ei gs and two-body contributions E 2gs . The former 
correspond to the kinetic and confinement energy parts. Remember first that the 1-body 
Hamiltonian is given by 



It is a simple exercise to obtain Ei gs by integration of expression (5) with this observable: 

Ei S s - g + 4 • 

More instructive is to prove that this equals 2 n r h rr , where h rr denotes the 1-body energy 
associated to each natural orbital. The calculation runs as follows: 

= (i + <A fv 2 r + 1)( ^ 



We have used (8) and the identities 

J f„(Q; P)P 2 dQdP = j f„(Q; P)Q 2 dQdP = r+ l - ] J2(2r + l)x r = 
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Now for the two-body contributions. The interelectronic repulsion potential in (2) is 



r\ 2 , so these contributions are of the form J2 rs c r c s L sr , with the L sr given by: 



fi 2 -co 2 



fsr(qi,pi)fsr(q2;p2)(qi - ^) 2 dq 1 dq 2 dpi dp 2 
h,(Qi)hr(Qi)(Qi - Q2) 2 h s {Q 2 )h r {Q 2 ) dQ 1 dQ 2 . 



Here h r are the usual harmonic oscillator eigenfunctions for unit frequency We consider the 
diagonal r = s first, whereby 



I2 2 -L0 2 

2-JUJJi 



1 

r+ 2 



and thus 



n r L r 



fi 2 - u 2 Ay/uJJI (uj + fJ,)(y/oJ + y/Ji ) 2 _ /i 2 - u 2 u + ji 



4/i 



2u 



We have used that the expected value of Q 2 for a harmonic oscillator eigenstate is r + \ 

and that the expected value of Q is zero. Notice that — — — < 1. Therefore, to fill up the 

2tu 

presumed energy gap (u 2 — fi 2 )/4:fi we have to "dig deeper". 

Now J h s (Q)h r (Q) dQ = for s ^ r, by orthogonality. A non-vanishing contribution of 
the off-diagonal part may then come (only) from the terms 



, ^-/i 2 . 



2y/UjJI 



1 2 



h r (Q)h r+1 (Q)Q dQ 



We compute: 



± 



2 2 00 

UT - /T 



/Ul/J, 



n r n 



r'tr+l 



r=0 



h r {Q)h r+l {Q)Q dQ 



±- 



/j 2 A^/u]i(yJuj — ,JJl) ( yfuj — y/JI\ r r + 1 



UjJL 



= ±(u 2 - ^ 



E 

r=0 



= ± 



J 2 — jJ 2 id — /J, 



4/i 2u 

Here we employ Yl^Lo( r + -O^ = (1 — x )~ 2 - The factor (r + l)/2 comes from the definition 
of the emission operators a) = (Q — iP)/y/2 (or the absorption operators), with a)h r = 
\Jr + 1 h r+1 . There is also an overall factor of 2 coming from two subdiagonals for each r. 

Obviously there are no other contributions. In order to minimize the energy we now have 
to choose minus signs whenever s — r ± 1, so our contention on the alternating sign rule in 
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the SLK functional for the harmonium model is proved; indeed, in this section we made no 
use of P gs whatsoever. The total energy comes out as 

u /i + u 2 / jJ, fi 2 — u 2 U) + n oj 2 — ji 2 oj — fl oj \i 
2 + 4 + 4/x 2oo 4> 2u~j~ ~ 2 + 2 ' 

as it ought to be. 

V. DISCUSSION 

That's all very well, the devil's advocate now concedes. But is it not rather baroque? At 
the heart of density functional theory there is the proof of existence of a functional yielding 
the ground state energy from d gs . We have managed to get it by a roundabout method 
equivalent to reconstructing the two-body state. Can't we just proceed directly? Yes, we 
can: the energy of the ground state is just (twice) the average energy of the Gibbs ensemble 19 
represented by (7). To wit, 

E gs = E[d gs ] = 2^(-±- [ + = t ^±E. 

Nearly all the exchange-correlation energy functionals currently used in quantum chem- 
istry trace their ancestry to that of Miiller. 9 ' 13-15 ' 20 Such approximations, written in our 
terms, are most often of the following form: 

1 °° f 
£ xc [d] = -- ^2 a (nj,n k ) / fj k {xi)V{\ri - r 2 \)f kj (x 2 ) dx x dx 2 , 

j,k=i J 

with integration both on spin and orbital variables. These are all actually recipes for d 2 . 
For the Midler functional a{n^n k ) = y/njuk. A handy list of alternatives is provided in 
Ref. 14. According to that reference, all of them (except for the Hartree-Fock functional) 
violate antisymmetry; nearly all of them violate the sum rule for d 2 ] as well as invariance 
under exchange of particles and holes for the correlation part. 

It is well known that the Midler functional is overbinding. Our own rigorous proof of 
this fact for real two-electron atoms 21 is much more transparent than the one in Ref. 13 and 
shows that definite positivity of the Coulomb potential does play a decisive role, whereas 
the extra minus signs in Midler's functional do not. For these very reasons the Midler func- 
tional's tendency to overcorrelate needs reexamination in harmonium. Differences between 
Coulomb and confining potentials are of course considerable; nevertheless, detailed analytic 
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comparison of the proposed functionals with the exact one remains an useful exercise, throw- 
ing some light, from our viewpoint, on the elusive correlation functional. This will be done 
elsewhere. 

Also, the remark at the beginning of this section pictures the harmonium "atom" as a 
system in equilibrium, with temperature depending on the strength of the interelectronic 
repulsion. Although matters are very different for confining potentials versus electrostatic 
ones, as well as for atoms with more than two electrons, it would seem to suggest that 
concentration on S xc [d] is a poor strategy. 

The Wigner function for the Hartree-Fock state for harmonium is given by the quasiprob- 
ability 

71" 

- J_ p -(R 2 +r 2 )y/(w*+n*)/2 -{P 2 +p 2 ) / ^ (^+^) /2 

so that the expressions (6) are replaced by the rather trivial 

A= (VWTW)j2 \ a _ 1= (i/^JTWm \ 

Coming to the correlation energy for harmonium: use of Phf gives 

x e - { RZ + r*)^^y-2 e - { P 2 + P 2 )/^T^)T2 d 3p d 3 R ^ ^ 

_ 3y/(a; 2 + /i 2 )/2 | 3 a; 2 3y/(a; 2 + /i 2 )/2 | 3/i 2 

4 ~~ 4y/{u 2 + fj?)/2 4 1 V ' U- 2 - ,, 2 )/2 

= 3vV + /i 2 )/2, 
and so the correlation energy is 

TP J? J? J" + V f^±Z\ 35 2 
E c :=E -E HF = 3[— 2 as ^0. 

March and coworkers have suggested a definition for the correlation energy density on con- 
figuration space on the basis of the Hartree-Fock wave function and the exact ground state 
for harmonium. However, relative momentum is as important as relative position in de- 
termining interelectronic correlation, and it seems more appealing to define a correlation 
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energy density on phase space, in the spirit of Rassolov 22 and of more recent work by Gill 
and coworkers. 23 We deal with this elsewhere. 
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